function [Subl]=substationLocate(Vlin,Ain,n)
%SUBSTATIONLOCATE A function used to locate substations
%	Input: substation numbers; Aij Connection Matrix; consumer locations with loads
%	Output:	substation locations (line: substation number  column:2)
%	August 7 2014, RIsoe Denmark, Qin Wang
%	Algrithm: Current idea to locate the substaion is: follow the linked power grid, at a fixed step length, find the largest points of it. 
if nargin ==3
	Vl=Vlin;
	Vltmp=Vl;
	A=Ain;
	Atmp=A;
	
	
	for m=1:n
        count=0;
		Vl=Vltmp;
		A=Atmp;
		sizetmp=size(Vl);		
	%----Get load/distance density ----
		radius = calculateRadius(Vl)/sqrt(n+1-m);
		for i=1:sizetmp(1)
	  		for j=i+1:sizetmp(1)
				if A(i,j)~=0 
					for k=0:10
						pointtmp = [Vl(i,1)+(Vl(j,1)-Vl(i,1))*k/10, Vl(i,2)+(Vl(j,2)-Vl(i,2))*k/10];
						density = 0;
						for l=1:sizetmp(1)
                            if (Vl(l,1)-pointtmp(1))^2+(Vl(l,2)-pointtmp(2))^2>0
							density = density + (sqrt((Vl(l,1)-pointtmp(1))^2+(Vl(l,2)-pointtmp(2))^2))/Vl(l,3);	
                            end
						end
						if count==0
							count=count+1;
							subltmp=[density k i j radius pointtmp(1) pointtmp(2)];
						else
							count=count+1;
							subltmp=[subltmp;[density k i j radius pointtmp(1) pointtmp(2)]];
                        end
					end
				end
			end
		end
		%----END get-----------------------%
		Subl(m,:)=subltmp(1,:);
		for i=2:count
			if subltmp(i,1)<Subl(m,1)
				Subl(m,:)=subltmp(i,:);	
			end
		end

		%----Delete inside points---%
		count=0;
		for i=1:sizetmp(1)
			if ((Vl(i,1)-Subl(m,6))^2+(Vl(i,2)-Subl(m,7))^2 <= radius*radius)
				Vltmp(i-count,:)=[];
				Atmp(i-count,:)=[];
				Atmp(:,i-count)=[];
				count=count+1;
			end
		end
		%----End Del----------------%
		
	end
else
	error('Error input numbers');
end
end
